Patterns and drivers of evapotranspiration in South American wetlands

Evapotranspiration (ET) is a key process linking surface and atmospheric energy budgets, yet its drivers and patterns across wetlandscapes are poorly understood worldwide. Here we assess the ET dynamics in 12 wetland complexes across South America, revealing major differences under temperate, tropical, and equatorial climates. While net radiation is a dominant driver of ET seasonality in most environments, flooding also contributes strongly to ET in tropical and equatorial wetlands, especially in meeting the evaporative demand. Moreover, significant water losses through wetlands and ET differences between wetlands and uplands occur in temperate, more water-limited environments and in highly flooded areas such as the Pantanal, where slow river flood propagation drives the ET dynamics. Finally, floodplain forests produce the greatest ET in all environments except the Amazon River floodplains, where upland forests sustain high rates year round. Our findings highlight the unique hydrological functioning and ecosystem services provided by wetlands on a continental scale.


Figure S3
| Spatial assessment of annual evapotranspiration across South American wetlands.Each boxplot refers to all annual ET (long term average; 2000-2015) pixels within a given land cover (wetland forest, wetland non-forest, open water, upland) and a flood category (most flooded pixels (M) and least flooded pixels (L), based on 50%-50% quantiles) for each wetland area.High importance of flooding on ET spatial pattern: wetlands where the most flooded pixels have higher ET (median values) than the least flooded pixels.Intermediate: the most flooded non-forest pixels have higher ET than the least flooded ones, but lower ET than forest or upland; or the most flooded forest pixels have higher ET than uplands.Little: no major difference among wetland and upland ET.S1).For the in situ measurements, monthly data were obtained by averaging all daily values in months with at least 75% of data available.In the map, towers located within wetlands are marked with a triangle, and those within uplands with a circle.For the BAN tower, SEBAL results are presented for two adjacent pixels (1 km far from each other), one located in a predominantly savanna and another in a mainly forest area, to show the differences in ET depicted by the model in this heterogeneous area.In this case, the tower likely receives contribution from both savanna and forest areas.

Supplementary note 1. Detailed description of the SEBAL algorithm
The SEBAL algorithm 8 estimates the instantaneous evapotranspiration (or latent heat) rate as the residual of the surface energy balance (Equation 1), using remote sensing and meteorological data (wind speed, specific humidity, surface air temperature and incoming shortwave radiation).
The main model premise is that the near-surface vertical air temperature difference is linearly related to the surface temperature 9 , and that there are two extreme pixels that characterize the landscape, namely the hot and cold pixels.At the hot pixel, the latent heat is assumed as zero so that all available energy (  − ) becomes sensible heat.Conversely, at the cold pixel all available energy becomes latent heat.
Our methodology adapts the SEBAL application by Laipelt et al 10 (originally using Landsat data) for MODIS imagery, and applied it within Google Earth Engine cloud computing environment (see SEBAL steps, as well as main input and output data in Figure S11).The instantaneous evapotranspiration rate estimated with Equation 1 is converted into 8-day evapotranspiration rate, which is the temporal resolution of the adopted MODIS products.This conversion is performed by assuming a constant evaporative fraction during the 8-day period.Finally, the 8-day values are averaged to yield monthly evapotranspiration.

Net radiation
Net radiation is calculated as: where  is the broad-band surface albedo,   the incoming short-wave radiation (  −2 ),   the incoming long-wave radiation (  −2 ), and   the outgoing long-wave radiation (  −2 ).
0 = 0.95 + (0.01) The broad-band surface albedo () is calculated following Tasumi et al. 12 : where   is a weighting coefficient and   the surface reflectance.
To calculate   , we used the equation suggested by ASCE-EWRI 13 : where  is the atmosphere pressure (),  the water in the atmosphere (),  ℎ the solar zenith angle over a horizontal surface, and   the unitless turbidity coefficient where   = 1 for clean air.
Atmospheric pressure is estimated as: where  is the elevation above sea level (m) obtained from the SRTM mission.
Water in the atmosphere is estimated according to Garrison & Adler 14 : where   is the actual vapor pressure () estimated as 15 : cos  ℎ and  2 equations are based on Duffie & Beckman 16 : where  ℎ is the solar zenith angle over a horizontal surface, and  the declination of the Earth,  = latitude of the pixel,  = hour angle and  = day of year.

Soil heat flux
Soil heat flux () is computed with the following equation, calibrated with remote sensing data and ground measurements at the flux towers: where   is the surface temperature (K).

Sensible heat flux
The following equation is used to estimate the sensible heat flux () 17 : where   is the air density (. −3 ),  the specific heat of air at constant pressure (. −1  −1 ) and  ℎ the aerodynamic resistance (s m -1 ) between two near-surface heights, 1 and 2, where 1 = 0.1 and 2 = 2 m above the zero-plane displacement height. is the temperature difference and represents a linear function of   , as proposed by Bastiaanssen et al. 9 : where  and  are empirically determined coefficients.
Since both  and  ℎ are unknown, SEBAL uses an iterative process.For the first iterative process,  ℎ is estimated assuming neutral stability 11 : ⁄ )  *  (17)   where  1 and  2 are the heights above the zero-plane displacement of the vegetation where  are defined,  * the friction velocity (m s -1 ) and  the von Karman's constant (0.41).
To estimate  * in the first iterative process, the following equation is used: where  200 is the wind speed (m s -1 ) at 200m and   the momentum roughness length (m).
200 is estimated as: where  * , is the friction velocity, and ℎℎ is assumed as 100 m 18 .
* , is estimated as: where   is the height of the GLDAS information and   the wind speed (m s -1 ).
is assumed as: where ℎ is the vegetation height (m), assumed as ℎ = 3 m.
In the iterative process,  ℎ is the near-surface temperature difference at the hot pixel and is calculated with the following equation: ℎ  ℎ ℎ  ℎ  (22)   where  ℎ ,  ℎ ℎ and  ℎ are the sensible heat, aerodynamic resistance and air density at the hot pixel, respectively.
For  ℎ and  ℎ the following equations are used: ℎ =   ℎ −  ℎ (24)   where   ℎ ,   ℎ and  ℎ are the land surface temperature, instantaneous net radiation and soil heat at the hot pixel, respectively.
To calculate the  and  coefficients of the linear relationship between   and , we consider   = 0 for the cold pixel (i.e.,  ℎ =0), which in combination with Equation 24 yields: The next steps are the final part of the first iterative process.The Monin-Obukhov length () defines the stability conditions of the atmosphere in the iterative process.This equation represents the height at which forces of buoyancy and mechanical mixing are equal 11 : where   is the air density (Kg m -3 ),  the specific heat o fair at constant pressure (J Kg -1 K -1 ),  * the friction velocity (m s -1 ),   the land surface temperature (K),  the von Karman's constant (0.41),  the gravitational acceleration (9.807 m s -2 ) and  the sensible heat flux (W m -2 ).
For the Equation 34, a value of 2 m is adopted instead of 200 m following the suggestion by Allen et al. 11 , in order to avoid numerical instability.
Finally,  ℎ and  * are estimated again using the values obtained from   and  ℎ : where   is momentum roughness length for each pixel, which was based on the following equation, calibrated with remote sensing data and ground measurements at the flux towers: The iterative process continues until the stability of stability conditions of the atmosphere is obtained and the absolute difference are lower than 0.1: With the end of iterative process, we have a stable value of  ℎ , and  is calculated using Equation 15.

Monthly evapotranspiration
The instantaneous latent heat is computed using the energy balance equation: The 8-day evapotranspiration ( 8− ) is then computed with the following steps.Firstly, the evaporative fraction (Ʌ) is calculated: Then, the latent heat of vaporization () (kJ kg -1 ) is estimated as: To further demonstrate the suitability of SEBAL in comparison to other ET algorithms, we performed additional analyses using GLEAM and MOD16 products, which are among the most widely used ET datasets (Figures S12 and S13).These two algorithms are not based on LST but mainly rely on vegetation indices to estimate land surface properties and feed ET algorithms.They are driven by meteorological reanalysis as an indicator of water availability (e.g., vapor pressure deficit and available radiation for MOD16, and precipitation and available radiation for GLEAM), relying on vegetation phenology as surface processes 27 .In addition, vegetation-based models are also dependent on biome properties look-up tables and land cover classifications, therefore most of the uncertainties reported are related to the low spatial resolution of meteorological reanalysis and land cover parameterization 28 , with higher uncertainties in areas with wetlands and high soil moisture 29,30 and higher accuracy over large and homogeneous areas 28 .Consequently, they are known to have difficulties to estimate differences between wetland and adjacent uplands, which can be depicted by LST-based algorithms such as SEBAL and other relatively similar algorithms like SSEBOP 10,23,31 .Using both MOD16 and GLEAM, we performed the same Budyko-like framework done in our analysis with SEBAL. Figure S13 clearly shows the lack of representation of the regional effect of inundated lands on ET for MOD16 and GLEAM as compared to SEBAL.Similarly, the comparison of the three different ET estimates with the measurements over the 10 flux towers (see Figure S13) shows that SEBAL also has relatively better performance metrics (R² and RMSE) than the other two datasets.These analyses support the choice of LST-like algorithms, as SEBAL, for estimating wetland ET at large domains.

Figure S1 |
Figure S1 | Climatology of average evapotranspiration in the wetland and adjacent upland regions.The value in each figure's title relates to the long term average ET for upland (Up) and wetland (Wet).Please note that ET data for Negro, Up and Down Amazon are not computed for the months from January to May because of persistent cloud cover which hampers the use of MODIS optical data.

Figure S2 |
Figure S2 | Long term flood fraction for South American wetlands.Each boxplot refers to long term (2000-2015) flood fraction values for all 25 km GIEMS-2 pixels within each wetland area.Red crosses are outlier values.

Figure S4 |
Figure S4 | Correlation among environmental drivers for the South American wetlands.Pearson correlation values are presented for the correlation among the following drivers for each wetland: flood fraction (FF), precipitation (P), leaf area index (LAI), net radiation (Rn), land surface temperature (LST), vapor pressure deficit (VPD) and wind speed (WS).*Significant correlations with P<0.01.

Figure S5 |
Figure S5 | Correlation between evapotranspiration (ET), ET to evaporative demand ratio (ET/E0), and environmental variables for South American wetlands.The assessed environmental variables are net radiation (Rn), leaf area index (LAI), flood fraction and precipitation (P).

Figure S6 |
Figure S6 | Climatology of evaporative demand, flood fraction, precipitation, and evapotranspiration to evaporative demand (ET/E0) ratio for the 12 wetlands.Please note that ET data for Negro, Up and Down Amazon are not computed for the months from January to May because of persistent cloud cover which hampers the use of MODIS optical data.

Figure S7 |
Figure S7 | Scatterplot between wetland-upland ET difference and flood fraction for each wetland.The Pearson correlation value is presented for each wetland.

Figure S8 |
Figure S8 | Climatology of wind speed (Ws), vapor pressure deficit (VPD), and evaporative demand (E0) for the 12 wetlands.E0 generally follows the VPD pattern, and Ws to a lesser extent.

Figure S9 |
Figure S9 | Validation of the SEBAL-based evapotranspiration (ET) climatology with in situ data from 10 in situ micrometeorological towers (three using the Bowen ratio and seven the eddy covariance method) located across South America.Because of lack of available MODIS data for many months in the tropical/equatorial regions, the SEBAL climatology is based on the whole time series (2000-2015), while the in situ data have different periods of availability (TableS1).For the in situ measurements, monthly data were obtained by averaging all daily values in months with at least 75% of data available.In the map, towers located within wetlands are marked with a triangle, and those within uplands with a circle.For the BAN tower, SEBAL results are presented for two adjacent pixels (1 km far from each other), one located in a predominantly savanna and another in a mainly forest area, to show the differences in ET depicted by the model in this heterogeneous area.In this case, the tower likely receives contribution from both savanna and forest areas.

Figure S10 |
Figure S10 | Validation of SEBAL-based monthly evapotranspiration (ET) with in situ data.All daily in situ measurements for the 10 flux towers presented in Figure S9 are plotted for each day of the year, and compared to the SEBAL-based 8-day estimates for the same period of data availability.
stable.Momentum and heat transport (  and  ℎ ) are computed using the following equations (following Paulson19 and Webb 20 ):For  < 0:

Figure S12 |
Figure S12 | Budyko-live assessment of various wetlands across South America, for three different remote sensing based ET datasets.